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1. Introduction 



Generally, the Hamiltonian of the n-component spin cubic model [TJ [2] can be written as 

H/k B T = -J2 [ K Si ■ Sj + M(s t ■ s 3 )% (1) 

<i,j> 

where T is the temperature, k B is the Boltzmann constant, and the sum on < i,j > includes 
all pairs of nearest-neighbor (NN) sites. The spin Si is an n-component vector located on the 
i-th site, namely Sj = (sn, s i2 , • ■ ■ , s in ), such that one and only one of the n components has a 
nonzero value ±1. This model is also called 'face-cubic model' because the spin can be regarded 
as a vector located at the center of an n-dimensional hypercube, and pointing to the center of 
one of the faces of the hypercube. 

This model combines the Potts degrees of freedom [31 H] with Ising degrees of freedom. 
Therefore the Hamiltonian can also be written as 

n/k B T = -J2 [KsiSj + M\5„ W1 (2) 

<i,j> 

where Sj = ±1 is an Ising spin, and cr, = 1, 2, ■ • • , n is a Potts spin. The partition sum of the 
cubic model can be written as 

Z = V* e ~ H/kBT = IT e (XSlSj+M)<5<T - CT j . (3) 

W,M {*}.{*} <m> 

For the special case coshi^ = e~ M , this partition sum can be mapped to (see Ref. [H [2] for 
details) 

Z = ^2(nx) N "n N % (4) 

{b} 

where x = eMs ' n ^ . This model is defined in terms of bond variables that can take the values 

n 

'absent' and 'present'. The bond configuration {b} is restricted to be Eulerian, which means 
that each site is connected to an even number of bonds. Nj, is the number of bonds, and iV c 
is the number of components. Typically, a component is a group of sites connected by bonds, 
but it can also be an isolated site. The model may thus be called 'Eulerian bond-cubic model', 
and one of its configurations is shown in Fig. [H Furthermore, n is no longer restricted to be 
an integer number in (j3J), it can be any real number. 

The Eulerian bond-cubic model has been studied by means of transfer matrix (TM) calcu- 
lations and a finite-size scaling analysis in Ref. [5]. In the region 1 < n < 2, critical points 



and three critical exponents were determined. To further investigate the nature of the phase 
transition and the critical behavior, especially the geometric properties of critical configurations 
of the model, one may also employ Monte Carlo simulations. However, the problem arises to 
design an efficient Monte Carlo algorithm for this model, in view of the nonlocal weight n Nc . 
A local update of the Metropolis type algorithm requires an nonlocal search to determine the 
change of the number of components. This, together with the critical slowing-down, would 
make the simulation very time-consuming in the critical region. 

For this reason, we make use of the 'coloring algorithm', which is proven to be useful for 
some models with nonlocal weights in their partition sums, such as the random-cluster |6] 
model and the O(n) loop model[7]. It was firstly proposed by Chayes and Machta[Hl [9], and 
was originally combined with the Swendsen-Wang algorithm [TO] to simulate the Potts model or 
the random-cluster model. Ding et al. [UJ extended the application of the 'coloring algorithm' 
to the simulations of the O(n) loop model on the honeycomb lattice using the Metropolis 
algorithm, and greatly improved the efficiency of the algorithm. Deng et al. [12] further 
proposed two efficient cluster algorithms by combining the 'coloring' trick with the Swendsen- 
Wang algorithm for loop models: the algorithm 1 and the algorithm 2. More applications of 
this 'coloring' trick can be found in Ref. [T3J EH]- 

By using the algorithm 2 [12J, the Eulerian bond-cubic model @ has been preliminarily 
simulated. However, only critical points were reported. In this work, we develop a variant 
of the algorithm 1 to simulate the Eulerian bond-cubic model on the square lattice. We pay 
attention not only to the thermodynamic properties of the model, e.g. the critical points, the 
thermal exponent and the magnetic exponent, but also to geometric properties such as the 
fractal dimensions of the critical configurations. 

The paper is arranged as follows: In Sec. [2we describe the cluster algorithm in details for 
the Eulerian bond-cubic model. In Sec. |3we describe the variables to be sampled and their 
finite-size scaling behaviors. Our numerical results are given in Sec. |4nd a summary is given 
in Sec. |5. 
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Figure 1: A typical configuration of the Eulerian bond-cubic model and the definition of 'dual 
sites', the sites i and j are the 'dual sites' of a and b, reversely, a and b are also the 'dual sites' 
of i and j. 

2. Algorithm 

The partition sum of the Eulerian bond-cubic model (jlj) can be written as 

2 = EE ( Nc ) (nx) N »(n - l) c "l c *, (5) 

{&} Cr =0 \ Cr I 

where c r is the number of 'red' components, and c g is the number of 'green' components, 
with c r + c g = n c . Here, each Eulerian bond configuration is decomposed into a number of 
'colored-Eulerian-bond' configurations. After this decomposition, the 'coloring algorithm' can 
be applied in the procedure of the Monte Carlo simulation of the model. We make use of the 
Ising spins sitting on the dual lattice to represent configurations: if two NN Ising spins Si and 
Sj on the dual lattice are different, the corresponding edge on the original lattice between S{ 
and Sj is occupied by a bond. The Eulerian bonds are precisely the domain wall of the Ising 
spins, there is a two-to-one correspondence between the Ising-spin configurations {S} and the 
Eulerian bond configurations {&}, see Fig. [TJfor example. In order to describe the algorithm 
more clearly, we define the 'dual sites': for an edge on the original lattice, the 'dual sites' of the 
two sites connected by the edge are the two sites on the dual lattice that sit at the two sides of 
the edge, reversely, the two sites connected by an edge on the original lattice are also the 'dual 
sites' of the two sites that sit at the two sides of the edge. See Fig. [TJ for example. Then the 
Swendsen-Wang type algorithm with 'coloring' trick for the Eulerian bond-cubic model can be 
described as: 



1. Start from an arbitrary Ising-spin configuration, which corresponds to an Eulerian bond 
configuration on the original lattice. 

2. Set color to the components of the Eulerian bond configuration. For each component on 
the original lattice, set it as green (active) with probability p — 1/n, or red (inactive) 
with probability 1 — p. 

3. Construct the Swendsen-Wang clusters by placing percolation bonds on the dual lattice. 

• For each pair of NN sites i and j on the dual lattice, a percolation bond is placed 
between them with probability p — 1 if not all of the colors of their dual sites are 
green. 

• If all of the colors of their dual sites are green and Si = Sj, connect them by a 
percolation bond with probability p = 1 — nx; otherwise, let the edge be vacant. 

• Each pair of NN sites on the dual lattice is considered to be in the same cluster 
if there is a percolation bond between them. These percolation clusters are called 
Swendsen-Wang clusters. 

4. Flip every Swendsen-Wang cluster with probability p = 1/2. 

5. Sample the variables of interest, erase the colors and restart at step 2. 

The algorithm can be modified to be Wolff- type [TB] by constructing only one cluster which is 
then flipped with probability p — 1 in step 3 and 4. 



3. Sampled variables and their finite-size scaling behavior 

A typical high-temperature (small x) configuration of the Eulerian bond-cubic model with 
a given n has only a few small components, as shown in Fig. [2](a). When x becomes larger, 
the typical configuration has more bonds and bigger components. A component that spans the 
whole lattice will emerge when x reaches or exceeds the critical point x c , see Fig. [2](b) and (c). 

The behavior of the components in these configurations is very similar to that of clusters 
in percolation phenomena [16] , so we call the spanning component a 'percolating component'. 
The percolation probability P s is defined as 

P. = (P) , (6) 



Figure 2: (Color online) Three typical configurations for n = 1.5 Eulerian bond-cubic model 
with system size L = 128: (a) at high temperature x = 0.3, (b) at critical point x = 0.444245, 
(c) at low temperature x = 0.6. Green means the spin on the dual lattice is '+', blue means 
'-'; the domain wall of the Ising clusters are exactly the Eulerian bond configurations. 



where P is 1 if there exists a percolating component on the configuration, otherwise. (• ■ ■) 
means the average over the canonical ensemble. For an infinite system, P s is 1 for x > x c and 
for x < x c , which is a G function. However, for a finite system, the value of P s changes 
continuously when shown in Fig. [3J 
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Figure 3: The percolation probability P s versus x for n = 1.5 Eulerian bond-cubic model with 
system size L = 256. 



On the other hand, the phase transition of this model can be described by in terms of the 
Ising spins on the dual lattice. When x is small, most Ising spins have the same sign, thus 
the system is in a long-range ordered state (ferromagnetic) and has a nonzero spontaneous 
magnetization, as shown in Fig. E(a). When x is large, the Ising spins will be in a disordered 



state (paramagnetic), and the magnetization will be zero, as shown in Fig. [2]Jc). The phase 
transition of the system is a ferromagnetic one. Concretely, the magnetization m is defined as 



rn 



with 



(M), 



12 S t 



(7) 



M 



V 



(8) 



where V = L d is the system volume and d is the dimension of the lattice. In the current paper, 
d = 2. Figure H] shows the magnetization m versus x for n = 1.5 Eulerian bond-cubic model 
with system size L = 256. Figure El and 0] give a general description of the phase transition 
of the model. 
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Figure 4: The magnetization m versus x for n = 1.5 Eulerian bond-cubic model with system 
size L = 256. 



The critical point can be determined by the percolation probability P s and the Binder ratio 
of m 

4 (M') ' { ' 

According to the renormalization theory, the percolation probability P s , the magnetization m 
and the Binder ratio Q display the following finite-size scaling behavior |17[ ITS] 

P s = Pj 0) + ai (x-x c )L yt +a 2 (x-x c ) 2 L 2yt + + + •••, (10) 

m = L ym - d [m + ai (x -x c )L yt + a 2 (x - x c ) 2 L 2yt + ... + b 1 L yi +b 2 L V2 + ■■•], (11) 
Q = Qo + ai{x - x c )L yt +a 2 (x-x c ) 2 L 2yt + ■ ■ ■ + hL yi + b 2 L y ' 2 + •••, (12) 



where y% > and y m > are the thermal exponent and magnetic exponent respectively, and 
yi, 2/2, • • ■ , < are the correction-to-scaling exponents. x c is the critical point and mo, dj, fej 
are unknown constants. We can see that P s and Q have a similar finite-size scaling behavior. 
Both (fTUl) and f[T2j) can be used to determine the critical point and the thermal exponent, as 
will be shown in more detail in Sec. |4. 
At the critical point x c , (TTTT) reduces to 



'in 



L ym - d {m + b x L m + b 2 L V2 + ■■■). (13) 



Both (ITT]) and (ITS]) can be used to determine the magnetic exponent y m , however, (1131) is used 
more in practice because it has fewer parameters thus is more convenient in a data analysis. 

Besides the critical exponents y t and y m , we are also interested in the fractal structure of the 
critical configurations. There are two fractals on the critical configuration: the largest Eulerian 
bond components and the largest Tsing cluster'. The Ising cluster is defined as a group of NN 
Ising spins in the same sign. We define the percolation strength P^ and the average size x° 
based on the components and Ising clusters. 

Pi - (14) 

i 

*" = £>■ (15) 

i 

where the superscript a = s, b or c. For a = s, Ci is the number of sites in the i-th component; 
for a = b, Ci is the number of bonds in the i-th component; and for a = c, q is the number of 
Ising spins in the i-th Ising cluster, is the size of the largest component or Ising cluster on the 
configuration. The subscript oo is used because only the largest component or the largest Ising 
cluster may have an nonzero fraction comparing to the system volume in the thermodynamic 
limit. The P^ can be considered as the order parameter of the phase transition, playing the 
role of the magnetization in the Ising model. The Greek letter x is used to denote the average 
size, because it has the property that is very similar to the magnetic susceptibility of the Ising 
model. 

P^ and x a have the finite-size scaling behaviors similar to that of the magnetization: 
P« = L y - d (P^ + ai (x-x c )L yt +a 2 (x-x c ) 2 L 2yt + --- + b 1 L y i +b 2 L y * + •■■), (16) 

x a = L 2y-d( Xo + a ^ x _ x j L y t +a ^ x _ Xc y L 2y t + ... + biLVl +b2L y 2 + ...^ ^ 



When a — s or b, y — yjj is the fractal dimension of the largest component of the critical 
Eulerian bond configuration, which can be viewed as the hull of the largest Ising cluster; when 
a = c, then y = y c is the fractal dimension of the largest critical Ising cluster. X°> a i an d 

hi are unknown constants, and the j/j < are irrelevant exponents. At the critical point x c , 
(USD and (HZD reduce to 

pa = £H(pM +6li » + ^ + ...) ) (1 8 ) 
X a = L 2 ^ d ( Xo + &i^ 1 +&2^ 2 + •••)• (19) 

We will use (TT8|) and (1TTJT) to determine y# and y c . 

4. Results 

Using the Swendsen-Wang type or Wolff type algorithm described in Sec. [2we do Monte 
Carlo simulations for the Eulerian bond-cubic model on the square lattice. We apply 10 4 
Swendsen-Wang/ Wo Iff cycles to equilibrate the system, and average over 2 x 10 7 samples, 
where each sample is taken after every 3 cycles. The sizes of the simulated systems range from 
L=8 to L=256. Figure |5] shows part of the data of P s versus x for n — 1.5 Eulerian bond- 
cubic model near the critical point x c . We fit the data according to ffTU]) using the nonlinear 
Levenberg-Marqurdt least-squares algorithm, which yields the thermal exponent y t = 0.747(6) 
and the critical point x c = 0.4443(2). 
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Figure 5: The percolation probability P s versus x for various system sizes of the n = 1.5 
Eulerian bond-cubic model. 




Figure |6] shows part of the data of Q versus x for n — 1.5 Eulerian bond-cubic model. 
Fitting the data according to f fT2|) . we obtain the thermal exponent yt = 0.748(3) and the 
critical point x c = 0.444245(8). We can see that results for y t and x c from the fit to the data 
for P s agree well with the ones from the fit to the data for Q. 
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Figure 6: The Binder ratio Q versus x for various system sizes of the n = 1.5 Eulerian bond- 
cubic model. 

According to Coulomb gas theory, for 1 < n < 2, the Eulerian bond-cubic model belongs to 
the same universality class of the critical O(n) loop model, with critical exponents [i~9 | l20 | I2T1 122] 
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(22) 


Vc 
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= 1 + -H , 

2 8g' 


(23) 



with n = — 2cos(7rg), 1 < g < 2. g is the Coulomb gas coupling constant. According to fl2TJ|) . 
we have yt = 0.748109 • • • for n = 1.5 Eulerian bond-cubic model, which is consistent with the 
numerical result. For the critical point of the Eulerian bond-cubic model on the square lattice, 
there is no exact result except for the cases n = 1 and n = 2. x c is ^2 - 1 for the case n = 1 
and 0.5 for the case n = 2 [2] - The critical points and some critical exponents of the Eulerian 
bond-cubic model are also determined by a finite-size scaling analysis based on numerical TM 
calculations in Ref. (2j, where it is found x c = 0.44424(1) and yt = 0.749(2) for the n = 1.5 
Eulerian bond-cubic model, which are in good agreement with our Monte Carlo results. 




We have also simulated the cases n = 1.0, 1.1, 1.25, 1.7, 1.75, and 2.0, both the Coulomb gas 
predictions for the critical O(n) branch and the numerical results for y t and x c are summarized 
in Tab. dJ 

Table 1: Thermal exponents and critical points of the Eulerian bond-cubic model on the square 
lattice. T = theoretical predictions of the critical O(n) branch, P s = Monte Carlo results from 
the fit to the data for P s , Q = Monte Carlo results from the fit to the data for Q, TM = results 
calculated by TM0. 



n 


yt(T) 


yt(Ps) 


y t (Q) 


2/t(TM) 


x c (P s ) 


Xc(Q) 


Xc(TM) 


<2o 


P(0) 
s 


1.0 


1.000000 


1.01(2) 


0.998(2) 


1.000000(1) 


0.41422(1) 


0.414214(2) 


0.4142135(1) 


0.8560(2) 


0.0651(3) 


1.1 


0.957313 


0.959(3) 


0.957(3) 


0.9572(5) 


0.41916(1) 


0.419155(4) 


0.419154(2) 


0.8359(2) 


0.0812(4) 


1.25 


0.887399 


0.890(4) 


0.888(2) 




0.42741(1) 


0.427404(3) 




0.8022(1) 


0.110(1) 


1.5 


0.748109 


0.747(6) 


0.748(3) 


0.749(2) 


0.4443(2) 


0.444245(8) 


0.44424(1) 


0.7454(3) 


0.1904(5) 


1.7 


0.600379 


0.598(4) 


0.604(6) 


0.595(5) 


0.4624(3) 


0.46213(2) 


0.46214(1) 


0.6686(3) 


0.294(2) 


1.75 


0.552482 


0.554(2) 


0.57(3) 




0.4679(4) 


0.46754(2) 




0.649(2) 


0.329(1) 


2.0 


0.000000 


0.52(3) 


0.48(3) 


0.50000(1) 


0.4998(3) 


0.50001(1) 


0.5000000(1) 


0.550(2) 


0.554(5) 



One should pay more attention to the case of n = 2. The thermal exponent y t is for 
the critical branch of 0(2) loop model, while it is 0.5 for n = 2 Eulerian bond-cubic model. 
For the case n = 2, the Eulerian bond-cubic model reduces to a special case of the Ashkin- 
Teller model[20l [23j [2U [25]. Our estimations of y t agree with the exact results y t = 0.5 for the 
Ashkin- Teller model. 

At the critical point, we sampled the magnetization m, percolation strength P^, P^, P^ 
and the average cluster sizes x c , X s \ X b - 

The log- log plot of the data for magnetization m versus system size L for n = 1.5 Eulerian 
bond-cubic model at the critical point is shown in Fig. [71 Fitting the data according to (fT3|) . 
we obtain the magnetic exponent y m = 1.7803(3), which is consistent with the theoretical 
prediction y m = 1.78054 • ■ -, given by f [2~T|) . 

The log-log plot of the data for P^ and \ c versus L are shown in Fig. [8] and Fig. [9] 
respectively. Fitting the data, we obtain the fractal dimension of the largest critical Ising 
cluster y c = 1.9195(10) (fit from P^) and y c = 1.9195(11) (fit from x c ), which are consistent 



with the theoretical prediction y c = 1.91989 ■ ■ •, given by fl23|) . 

Fitting the data of P^, x s , P^, and we obtain the fractal dimension of the largest critical 
component y H = 1.410(5) (fit from y H = 1.405(4) (fit from x s ), Vh = 1.410(4) (fit from 

P^), and yn = 1.41(1) (fit from x b )- All of them are consistent with the theoretical prediction 
y H = 1.4064 • • given by (j22j|. 

All the numerical results and the theoretical predictions of y m , y c and yn are listed in Tab. 
[2] and Tab. [3j From these tables, we can see that the fractal dimension y c of the largest critical 
Ising cluster of the n = 2 Eulerian bond-cubic model is the same as the critical 0(2) value. 
However, the fractal dimension of the largest critical Eulerian bond component is yu = 1.625, 
which is obviously different from the critical 0(2) value yjj = 1-5, given by (122]) . Also the 
magnetic exponent y m is different from the critical 0(2) value 1.5. These results agree well 
with the exact results y m = 1.625 = yu for the Ashkin- Teller model [221 EH and again 
show the difference between the cubic symmetry and the O(n) symmetry in the case n = 2, 
when the cubic anisotropy becomes marginally relevant. 

Table 2: The magnetic exponents y m and fractal dimension y c of the Eulerian bond-cubic model 
on the square lattice. T = theoretical predictions of the critical O(n) branch, MC = Monte 
Carlo results, = Monte Carlo results from the fit to the data for P^, \ c = Monte Carlo 
results from the fit to the data for x c , TM = results calculated by TM[2]. 



n 


Z/m(T) 


y m (MC) 


Vm(TM) 


Ve(T) 




ydx c ) 


1.0 


1.875 


1.8751(2) 


1.87501(1) 


1.94792 


1.9476(3) 


1.9477(3) 


1.1 


1.85899 


1.8590(1) 


1.85895(5) 


1.94257 


1.9422(3) 


1.9423(3) 


1.25 


1.83277 


1.832(1) 




1.93436 


1.933(2) 


1.933(2) 


1.5 


1.78054 


1.7803(3) 


1.7805(5) 


1.91989 


1.9195(10) 


1.9195(11) 


1.7 


1.72514 


1.724(1) 


1.726(2) 


1.90702 


1.906(1) 


1.909(3) 


1.75 


1.70786 


1.709(2) 




1.90347 


1.902(2) 


1.90(1) 


2.0 


1.5 


1.6249(4) 


1.62500(1) 


1.875 


1.875(1) 


1.8750(1) 



Table 3: The fractal dimension yn of the Eulerian bond-cubic model on the square lattice. T 
= theoretical predictions of critical O(n) branch, P^ = Monte Carlo results from the fit to the 
data for P^, \ s — Monte Carlo results from the fit to the data for \ s , an d so forth. 



n 


Vh(T) 


Vh{PL) 


Vh(x s ) 


y H (P b J 


VH(x b ) 


1.0 


1.375 


1.376(1) 


1.373(4) 


1.378(4) 


1.372(4) 


1.1 


1.3803 


1.3806(9) 


1.380(5) 


1.382(4) 


1.379(3) 


1.25 


1.3891 


1.391(2) 


1.390(2) 


1.393(4) 


1.389(1) 


1.5 


1.4064 


1.410(5) 


1.405(4) 


1.410(4) 


1.41(1) 


1.7 


1.4249 


1.426(4) 


1.43(2) 


1.42(1) 


1.430(6) 


1.75 


1.4307 


1.43(1) 


1.431(12) 


1.427(11) 


1.43(1) 


2.0 


1.5 


1.625(1) 


1.6253(5) 


1.625(1) 


1.624(1) 




Figure 7: The log- log plot of the magnetization m versus system size L for n = 1.5 Eulerian 
bond-cubic model. The error bars are much smaller than the size of the data points. 
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L 

Figure 8: The log- log plot of the percolation strength versus system size L for n — 1.5 
Eulerian bond-cubic model. The error bars are much smaller than the size of the data points. 




Figure 9: The log-log plot of the average size x c versus system size L for n = 1.5 Eulerian 
bond-cubic model. The error bars are much smaller than the size of the data points. 



5. Summary 



We simulated the Eulerian bond-cubic model on the square lattice using an efficient cluster 
algorithm. Two fractal dimensions of the critical configurations as well as the critical points, 
the thermal and magnetic exponents of the model are determined by means of a finite-size 
scaling analysis. The two fractal dimensions are for the first time obtained numerically. The 
estimations of the critical points and the thermal and magnetic exponents are in good agreement 
with those obtained by means of TM calculations [2] for several values of n, including the case 
n = 2. Our results for all critical exponents are consistent with the Coulomb gas predictions 
of the critical O(n) branch for n < 2. But, for n = 2, the thermal exponent, the magnetic 
exponent and the fractal dimensions yu are different from the critical 0(2) values, and the 
model reduces to a special case of the Ashkin- Teller model. Our study confirms that the phase 
transition of the Eulerian bond-cubic model belongs to the critical 0(n) universality class for 
n < 2. The cubic anisotropy is irrelevant for n < 2, but becomes marginal when n — 2. 
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